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Abstract 

This work describes two algorithms for computing the angular rate and attitude in case of 
a gyro and a Star Tracker failure in the Microwave Anisotropy Probe (MAP) satellite, which was 
placed in the L2 parking point from where it collects data to determine the origin of the universe. 
The nature of the problem is described, two algorithms are suggested, an observability study is 
carried out and real MAP data are used to determine the merit of the algorithms. It is shown that 
one of the algorithms yields a good estimate of the rates but not of the attitude whereas the other 
algorithm yields a good estimate of the rate as well as two of the three attitude angles. The 
estimation of the third angle depends on the initial state estimate. There is a contradiction between 
this result and the outcome of the observability analysis. An explanation of this contradiction is 
given in the paper. 

Although this work treats a particular spacecraft, the conclusions have a far reaching 
consequence. 


Introduction 

The Microwave Anisotropy Probe (MAP) Satellite was launched at 15:46:46 EDT on June 
30, 2001 aboard a Delta 11-7425- 10 launch vehicle. After three phasing loops MAP passed by a 
swingby point where it was injected by the moon’s gravity towards the L2 point. On October 1, 
2001, following a three month journey, MAP arrived safely at its permanent observing station 
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near the L2 Lagrange Point, 1.5 million km from Earth, a quasi-stable position a million miles 
from Earth in the direction opposite the Sun. 

MAP used two two-axis gyroscopes to measure its angular rate vector. One input axis of 
the first gyroscope was aligned along the body x-axis and the other input axis was aligned along 
the body z-axis. The input axes of the second gyro were aligned along the body y and z-axes. 
Therefore, if one gyro failed then the rate about two body axes is still measured; namely, about 
the z-axis and one of the other two. MAP also had two Digital Sun Sensors that measured one 
vector; namely, the direction to the Sun. In addition, the spacecraft also carried two Star 
Trackers. However, when MAP performed the phasing loops it passed through the Van- Allen 
radiation belts, where the star trackers were inoperative. 

The problem handled in this paper is as follows. Suppose that indeed one gyro failed and 
the star trackers were inactive when MAP passed through the Van-Alien belts. Was it still 
possible to reconstruct the full angular rate and the orientation of MAP with the one operating 
gyro and a measurement of the Sun direction? It should be noted that although this question 
concerns a particular mission, its implication is more general. 

We addressed this problem by defining two possible filtering models and performing a 
nonlinear observability analysis. Then we performed simulations of the spacecraft motion and ran 
the suitable filters to verify the results of the analysis. In the next section we describe the angular 
rate at which the problem is examined. This section is followed by a description of the examined 
filters. In the section that follows we present the observability analysis that was carried out. In the 
next section we present simulation test results and then, in the last section, we present our 
conclusions. 
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Body Angular Kinematics 


During the critical part of the trajectory; that is, when MAP passed through the radiation 
belts the spacecraft was kept at an inertial hold where the angular rate was nominally zero. In this 
mode the spacecraft attitude control system acted as a regulator where the gyros nominally read 
no input. In order to examine our suggested filters under more stringent conditions we chose to 
examine the filters under non-zero nominal angular rate profile. As a baseline rate profile we 
chose the MAP science collection rate profile at L2. 

At the L2 observation station the angular motion of MAP consists of two rotations as 
described in Fig. 1. One rotation, at 360 degrees an hour, is about the anti-Sun line, and the other, 


Sun 

X direction 

0.464 rpm 

* 1 .00 rph 



Fig. 1 : The Geometry of the Angular Motion Components of the MAP Spacecraft 
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at 0.464 RPM, is about the satellite z-axis. The body z-axis is tilted by 22.5 degrees with respect 
to the anti-Sun line. The body angular rate and attitude definitions are presented in Fig. 2 in terms 


Sun 

f] 



Fig. 2: Coordinate Systems Involved in the Description of the MAP Total Angular Velocity 

of three Euler angles. In order to describe the spacecraft motion and the associated measurements 
we define four coordinate systems. The first system is the Sun coordinate system (designated by 
sun). The center of this system is at the origin of the spacecraft coordinates, its z-axis points at the 
Sun direction and its x-axis, which is turning, is assumed to point at some arbitrary direction in 
space, which we define as zero y. A rotation by the angle -y about the z-axis yields a first 
interim coordinate system, which we denote by /. A second rotation by the angle 6 about the y- 
axis of this interim coordinate system yields the second interim coordinates, which we denote by 
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2. Finally, a rotation about the z-axis of the latter system results in the body coordinates, which we 
denote by b. For sake of simplicity, only the body system z-axis is shown in Fig. 2. 






Fig. 3: The Piogram Describing the Geometry and Rates of the MAP Spacecraft 


MAP rotates about two axes. One rotation, at the rate of - y , is about the Sun line, which 
we defined as the Z sun axis, and the other, at the rate of f) , is about the body z-axis. In order to 
express the sum of these rotation rates in body coordinates we use the Piogram [1] of Fig. 3. From 
the Piogram we obtain 


co x = ysin6-cosrj 
co y = -ysin8 -sinri 
co z = -ycos6 + f] 


and the transformation matrix from the Sun to body coordinates is: 


D 


sun 

b 


cycbcr) + sysrj - sycbcr) + cysrj - s8cr| 
- cycSsr) + sycr) syc5sr) + cycrj s8st] 

cysS - sys8 c8 


0 ) 


( 2 ) 


The corresponding quaternion can be extracted best when using the newly derived algorithm of 
[2]. The angular rates and angles are as follows: 
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y = 1 rev/hour = 1 • 2n/3600 rad/sec 

and 

y = yt + 0 rad = yt rad 

(3-a) 

f) = 0.464 rpm = 0.464 • 27i/60 rad/sec 

and 

q = qt + 0 rad = qt rad 

(3-b) 

8 = 22.5° =22.5 -ti/ 180 rad 



(3.c) 


Upon examination of the MAP trajectory, the Sun direction and thus the Sun coordinates can be 
considered stationary in inertial space relative to the rotation of MAP with respect to the Sun 
coordinates. Consequently, we may consider the rate components of Eq. (1) as inertial rates 
resolved in body coordinates. 

The Sun Sensor Measurements 

We wish to know the relationship between the sun sensor measurements and the angles y , 
8 and r| that determine MAP’s attitude. To do that we turn to Fig. 4 which depicts the 
orientation of the sun sensor boresight in body coordinates. The sun sensor reads the tangent of 
the angles a and (3 . We wish to express the readings of these tangents as functions of the angles 
y, 8 and q defined in Fig. 2. 



Fig. 4: Definition of the Sun Direction 
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To do that we define (Z sun ) b and (Z sun ) sun . The column matrix (Z sun ) b contains the components 
of the vector Z sun when resolved in the body coordinates, and, similarly, the column matrix 
(Z sun ) sun contains the components of the vector Z sun when resolved in the Sun coordinates. We 
denote the components of the former by a x , a y and a z . It is easy to see that the components of 


(Z su „ ) sun are 0,0, and 1. Thus 


and 

Obviously 


(Z sun )£=[a x a y a z ] 

( Z sun)L=[° 0 1] 

(Zsun )b= Dr (Z sm ) sun 


(4.a) 

(4.b) 

(5) 


From the last three equations and Eq. (2) we obtain 


V 


-s5cq 

'O' 

a y 

= 

x x s5sq 

0 

. a z_ 


c5 

! 


hence 


a x = -s8cr) 


a y = s5sr) 


therefore 


a z = c5 


( 6 ) 


(7) 


a 

tan p = — = - tan 5 • cr| 
a z 

( 8 ) 

a v 

tana = — = tano-sq 

a z 

Note that tana and tanp are independent of y. An inspection of Fig. 2 reveals that this is an 
expected result. 

Since tana and tanp are measured quantities, the last two equations contain only two 
unknowns, therefore we can solve for the angles 5 and q. For this reason even if the 
measurements contain unbiased errors we can still find these two angles by proper filtration. In 
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order to fully know the spacecraft attitude we need to know y too. If the gyros mounted on the 
Sun coordinates z-axis introduce no bias, then a proper integration of their measured outputs will 
yield y and thus we can obtain the spacecraft attitude from the measurements obtained from 
these gyros and the sun sensor. If the gyro measurements contain a bias then we ask ourselves 
whether the system is observable; that is, whether the spacecraft angular dynamics (Euler) 
equation may add information with which we can still observe y and thus completely determine 
the attitude. To answer this question we have to design an estimator and examine the 
observability of the system used by the estimator. This is done next. 


Filter I 


The Filters 


Dynamics Model 

In this estimator (filter) we convert the sun sensors measurements to a vector 
measurement, which is then connected to the quaternion of the spacecraft. The dynamics model of 
Filter I was developed in [3]. Using 


q 4 

-q 3 

q 2 

* 

q 4 

-q. 

-q 2 

q. 

q 4 

_-qi 

-q 2 

-q 3 


define the system dynamics as; 


CO 


A 



I"[(Ico+h)x] 0 

±Q 


o' 

(0 

+ 

'r'(T-h)' 

+ 


0 _ 

_q_ 


0 


A. 


(9) 


( 10 ) 


where [gx] is the cross product matrix of a general vector g, w M accounts for the inaccuracies in 


the modeling of the spacecraft angular dynamics, and w q accounts for modeling errors of the 
quaternion dynamics. 
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Measurement Equations 

The measurements are the tangent of the angles a and j8 , measured by the sun sensor and shown 
in Fig. 4, and the z and either the x or the y component of the body rate measured by one of the 
spacecraft gyros. The measured quantities are contaminated by noise so they are not the reference 
values but rather, a m , P m , co z m and either co x m or co ym , depending which gyro failed. 


Sun Sensor Measurement 

Using the two measured angles we compute a unit vector in the assumed sun direction as 
follows [4, p. 226] 


b 


m 


l 

V tanaf + tanp^T 


tan Pm 

tana m 

1 


( 11 ) 


In reality tanp m and tana m contain noise; however in the filter we approximate the b m 
measured vector as the true vector to the Sun and a simple additive zero mean white noise; that is, 

b„=b + v b (12) 

where b is the nominal unit vector in the true direction to the sun. The logic behind this 
simplification is that if in tests the filter yields satisfying results then the approximation is 
justified. Let r denote the measured Sun-vector expressed in the reference coordinate system 
(this vector is taken from the almanac). The relationship between the two vectors is expressed by 

b = Dr (13) 

where D is the direction cosine matrix that transforms vectors from the reference to the body 
coordinates. From the last two equations we obtain 

b m =Dr + v b (14) 

It is well known that D is the following function of the quaternion elements [4, p. 414] 
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(15) 


D = 


^ 2 2 2 , 2 

^2 Qj +< l4 

2(q,q 2 -q 3 q 4 ) 
2(q x q 3 + q 2 q 4 ) 


2(q,q 2 +q 3 q 4 ) 
-q?+q 2 -q 3 +q 4 
2(q 2 q 3 — qiq 4 ) 


2(q 2 q 3 q 2 q 4 ) 
2(q 2 q 3 +q,q 4 ) 
-qf-q 2 +q 3 +q 4 


when the last expression is substituted into Eq. (14) we obtain 


K = 


(q, 2 -q 2 -q 3 +q 4 )q + 2(q,q 2 +q 3 q 4 )r 2 +2(q,q 3 -q 2 q 4 )r 3 
2(q,q 2 -q 3 q 4 )r, +(-q, 2 +q 2 -q 3 +q 4 )r 2 + 2(q 2 q 3 +q,q 4 )r 3 
2(q,q 3 +q 2 q 4 ) r i +2(q 2 q 3 -q,q 4 )r 2 +(-qf -q 2 +q 3 +q 4 )r 3 


+ v fc 


(16) 


or 


b m = 


qf r i +q,q2 r 2 +q.q 3 r 3 -q 2 2r i +q,q 2 r 2 -q 2 q 4 r 3 -qft +q 3 q 4 r 2 +qiq 3 r 3 +q 2 r, +q 3 q 4 r 2 -q 2 q 4 r 3 
qiq 2 h -qf r 2 +q.q 4 r 3 +q.q 2 r, +q 2 r 2 +q 2 q 3 r 3 ~q 3 q 4 r, -q 3 r 2 +q 2 q 3 r 3 -q 3 q 4 r, + q 4 r 2 +q,q 4 r 3 
q,q 3 r, -W: -qfr 3 + q 2 q 4 r, +q 2 q 3 r 2 -q 2 r 3 +q,q 3 r, +q 2 q 3 r 2 +q 3 r 3 +q 2 q 4 r, -q,q 4 r 2 +q 4 r 3 


+ v. 


which can be written as 


..(17) 



" qi r,+q 2 r 2 +q 3 r 3 

-q 2 r.+qir2-q 4 r 3 

-q 3 r i+q 4 r 2+q. r 3 

q 4 ri+q 3 r 2 -q2 r 3 ' 

V 

q 2 

q 3 

_q 4 _ 

b m = 

q 2 r,-qir 2 +q 4 r 3 

q,r 1 +q 2 r 2 +q 3 r 3 

-q 4 q-q 3 r 2 +q2 r 3 

-q 3 ii+q 4 r 2 +q,r 3 


_q 3 r 3 -q 4 r 2 -q t r 3 

q 4 r,+q 3 r 2 -q 2 r 3 

q.q+q^+q^ 

q 2 r,-q,r 2 +q 4 r 3 _ 


+ Vi 


Let 


H b (r,q) = 


q.r.+q^+q^ 

-q 2 ri+q.r 2 -q 4 r 3 

-q 3 r,+q 4 r 2 +q,r 3 

q 4 r i+q 3 r 2-q 2 r 3 ’ 

...(18) 

q2q-qi r 2+q 4 r 3 

q,ri+q 2 r 2 +q 3 r 3 

-q 4 q-q 3 r 2 +q 2 r 3 

-q 3 ri+q 4 r 2 +q,r 3 

(19) 

q 3 r .-q 4 r 2 -q.r 3 

q 4 r .+q 3 r 2 -q 2 r 3 

qi r i +q 2 r 2 +q 3 r 3 

q 2 r i-qi r 2+q 4 r 3 . 



then Eq. (18) can be written as 
or 


*> m =H b (r,q)q + v b 


b m =[0 H b (r,q)] 


+ v b 


( 20 ) 

( 21 ) 


This is the measurement equation associated with the sun sensor measurement. 

Gyro Measurement 

Let us assume that the gyro which measures rates along the x and z axes continues to 
operate. The x and z axes rate measurements can be written as 
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“x.m = + V cox 

G>z,m =®z +V , Z 


( 22 ) 


where co x and co z are the true rates, and v mx and v (0Z are scalar zero-mean white-noise 


measurements. Define 


H,„ = 


1 0 0 
0 0 1 


then from the last two equations we obtain 


(23) 


co r 


CO. 


z,m 


=[«. <U 


V 

L coz J 


(24) 


We denote the measurement vector by u m where, clearly, co* = [co xm co zm ]. Similarly we 


denote the noise vector as by v m where = [v^ v wz ] . Obviously, if the other gyro is the one 

that is still operating, then the x component in Eqs. (22 - 24) has to be replaced by the y. 

The Combined measurement 

Combining Eqs. (21) and (24) we obtain 



0 


3 x 3 


ff*(r.q) 


0 


2 x 4 J 


CO 

q 



(25) 


This is the combined measurement model for simultaneous measurement of the operating two- 
axes gyro and the sun sensor. 

Filter II 

Dynamics Model 

In this filter we make use of the kinematics rule governing the rate of change of the Sun 
vector. It is known that 


Dr = b + w x b (26.a) 

hence 

b = [bx]w + Dr (26.b) 


11 



As mentioned earlier, the vector r is known from the almanac, therefore r is known too. 
Consequently we can treat Dr as a deterministic forcing function. In our case r is the vector from 
the spacecraft to the Sun. For the MAP trajectory this vector changes very slowly. Therefore we 
can neglect the term in r and use 


b = [bx]co (26.c) 

Augmenting Eq. (26.c) with the dynamics of Filter I yields 


d> 


'r'[(Iw +h)x] 

0 

o' 

u 


'r'(T-h)" 



q 

= 

IQ 

0 

0 

q 

+ 

0 

+ 

W q 

b 


[bx] 

0 

0 

b 


Dr 


- W b_ 


(27) 


The state vector in Eq. (27) includes the Sun vector b even though b is a measurement. Indeed, 
in order to formally compute the rate and attitude one can numerically differentiate the measured 
b x and b y , and use them in the following equations (which are the components of Eq. (26.c)) 


b x =co y b x -co x b y 
b y = © z b x -co x b z 


(28) 


together with the gyro-measured co z , to solve for the two unknowns co x and co y . Then use these 
co x and co y and the measured co z as ‘measurements’ in an estimator whose dynamics model 

consists of the first two matrix equations in Eq. (27). This, however, required numerical 
differentiation, which introduces noise. To avoid differentiation we apply the estimation approach 
[5]. This can be accomplished by adding Eq. (26.c) to Eq. (10), the differential equations of the 
variables that we want to estimate, thus forming Eq. (27). At this point a Kalman filter can be 
utilized to estimate those variables. 


12 


Measurement Equations 

The measurement equations of Filter II are quite simple. The vector measurement of Eq. 
(21) is replaced by the simple equation 


b m =[0 0 I,] 


CO 

q 

b 


+ V„ 


(29) 


Eq.(24) stays the same; therefore, the augmented measurement equation for Filter II becomes 

co 


L W mJ 


0 0 I, 

H. 0 0 


v 

L%J 


(30) 


Observability Tests 

Before carrying on an observability test it is necessary to specify the kind of filter that 
was used in this analysis. The filter that was used here is the PSELIKA (Pseudo-Linear Kalman) 
Filter [6], The idea on which this filter is based is as follows. Suppose that the non-linear 
dynamics and measurement equations of a non-linear system can be decomposed in the following 
way 

x = F(x)x + Bu + w (31. a) 

y = H(x)x + v (31 .b) 

If F(x) is close enough to F(x) and H(x) is close enough to H(x) (note that we do not assume 

that x is necessarily close to x ) then instead of the unknown matrices F(x) and H(x) we can 

use respectively F(x) and H(x) which are known along the estimation process and apply the 

linear Kalman filter algorithm to the measurements using the models 

x = F(x)x + Bu + w (31. c) 

y = H(x)x + v (31.d) 

The observability analysis requires F(x) andH(x). The x for F(x) and H(x) was the simulation 

truth model state. The logic for this choice of x is as follows. If when using F(x) and H(x) in the 
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analysis the latter indicates that the system is not observable under these favorable conditions 
then there is no sense in trying to estimate the state vector x of the system. If, on the other hand, 
the analysis indicates that the system is observable, then although observability and therefore 
estimability are not assured, it makes sense to design an estimator and attempt to estimate x . 

The observability of the system models used by the two filters was examined using the 
following analysis. Consider O k , the transition matrix, which corresponds to F, the dynamics 
matrix, of either Eq. (10) or Eq. (27). This matrix transforms the system state vector at time t k , 
x k , to x k+1 , the state at time t k+1 . If at a certain time-point, t m , the initial state, denoted by x 0 , 

can be computed, then, for our purposes, the system is observable. Adopting the common 
approach to proving complete observability of a discrete linear system, we express the first m 
measurements as follows 


yj = H j( x j) 


n®^.) 

i=0 


x 0 j = 0, 1 m 


(32) 


where O 0 (x 0 ) = I . Form the matrix equation 


’y 0 " 


H 0 (x 0 ) 

y. 



y 2 

- 

H 2 (x 2 )d> 1 (x 1 )d> 0 (x 0 ) 

_y m _ 


.H m (x m )O m . 1 (x m . 1 )---O 0 (x 0 )_ 


Note that in Filter II H is not a function of the state (see equations 23 and 30). That is, the 
measurement equation is linear from the start. If there are n independent rows in the right hand 
side matrix (the observability matrix) in Eq. (32) where n is also the number of states, then x 0 can 
be solved for and the state is observable. 


14 


The observability test of Filter I reveals that the rank of the observability matrix for this 
filter becomes 7 (see Fig. 5), and since the filter state is 7 too, the system is observable. It means 
that there is a good chance that Filter I can estimate both the missing angular rate component as 
well as the attitude of the spacecraft when one gyro fails. We note though that the measurement 
equation of Filter I is only pseudo-linear and not truly linear. This may cause a problem even 
though the observability test indicates that the state is completely observable. 

Filter I 



Fig. 5: Rank Development of the Observability Matrix of Filter I. 


Filter II 



Fig. 6: Rank Development of the Observability Matrix of Filter II. 
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The observability matrix of Filter II becomes 6 and stay at 6 (see Fig. 6) whereas the size 
of the state vector is 9. This indicates that the system is not completely observable. It means that 
we cannot fully estimate the missing rate component and the attitude. However, maybe their 
estimated values are close enough to the truth. Therefore both filters had to be tested. Simulation 
results of the two filters are presented in the next section. 

Simulation Test Results 

Simulations were run for different initial conditions and noise levels. Plots that resulted 
from a typical run of Filter I are presented in Figs. 7 - 12. For this run the initial y angle was 60 
degrees, 6 was 35 degrees and t| was 0. The initial estimates were zero for all three angles. The true 
Euler-angle rates were 2-(2 ti)/ 3600 rad/sec for y, zero for 5 and 0.928 • (2 tt)/ 60 rad/sec for f|. 
According to Eq. (1), using these initial rates and angles, the true initial angular rate components 
were as follows: to =0.002 rad/sec ct>„ =0 rad/sec and ©_ =0.094 rad/sec. The initial estimated 

values of 



Time [sec] 

Fig. 7: The ^-component of the angular velocity, its estimate, and its estimation error for Filter I. 
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the angular velocity were zero for all three components. The one-sigma gyro measurement noise was 
0.01 deg/sec and that of the sun sensor was 0.25 degrees. In each figure we see the true state 
superimposed on the estimate of the state. Underneath this plot there is a plot of the estimation error. 
In this run y is 




Fig. 8: The y-component of the angular velocity, its estimate, and its estimation error for Filter I. 



Time [sec] 


Fig. 9: The z-component of the angular velocity, its estimate, and its estimation error for Filter I. 
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Time [sec] 

Fig. 10: The y angle, its estimate, and the estimation error for Filter I. 



0 20 40 60 80 100 120 140 160 180 200 


Time [sec] 

Fig. 11: The 8 angle, its estimate, and the estimation error for Filter I. 
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Fig. 12: The r| angle, its estimate, and the estimation error for Filter I. 

estimated quite well, but this is not always the case. The observability of this angle is quite poor. 
Filter II was also run but, as expected from the observability test, the three attitude angles were 
not observable. The estimates of the angular rate components resemble those obtained using Filter 
I but the estimates of the three angles were quite poor as seen in Figs. 13-15. 




Time [sec] 


Fig. 13: The y angle, its estimate, and the estimation error for Filter II. 
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Time [sec] 


Fig. 14: The 5 angle, its estimate, and the estimation error for Filter II. 




Fig. 15: The r| angle, its estimate, and the estimation error for Filter II. 


Conclusions 

In this paper we considered a spacecraft that uses two two-degree of freedom gyros and 
two vector measuring devices for measuring its rate and its attitude. We investigated the 
possibility of estimating the spacecraft full rate vector and attitude if one gyro and one vector- 
measuring device fail simultaneously. In our analysis we proposed two filters for estimating the 
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desired variables and carried out an observability analysis for both filters. The observability study 
indicated that the variables were completely observable when the filter designated as Filter I was 
used, and was not completely observable when the other filter, designated as Filter II, was used. 
We tested our conclusions through simulations and found that, as expected, Filter II was unable to 
estimate the attitude angles. However, it was also found that, contrary to the results of the 
observability analysis, when using Filter I the estimation quality of one of the attitude angles 
depended on the initial errors. This is explained by the fact that the observability matrix was state 
dependent. Therefore only when we ran the test using the correct state; namely, the nominal 
angular rate and Euler angles, was the system found to be observable. In reality, however, the 
existing state vector was the estimated state vector. Therefore as long as the initial estimate was 
close enough to the true state, the simulation conditions matched those of the analysis and the 
observability test results predicted well the performance of the filter. On the other hand, when the 
initial estimate of the state vector was far from the true vector the analysis was not valid anymore. 
The simulated data used in the analysis were those of the MAP spacecraft attitude and rate profile 
at the L2 Lagrange point. 
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